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ABSTRACT - The Global Positioning System (GPS) has become a standard 
method for low cost onboard satellite orbit determination. The use of a GPS 
receiver as an attitude and rate sensor has also been developed in the recent past. 
Additionally, focus has been given to attitude and orbit estimation using the 
magnetometer, a low cost, reliable sensor. Combining measurements from both 
GPS and a magnetometer can provide a robust navigation system that takes 
advantage of the estimation qualities of both measurements. Ultimately a low 
cost, accurate navigation system can result, potentially eliminating the need for 
more costly sensors, including gyroscopes. 

1 - INTRODUCTION 

This work presents the development and initial testing of a unified navigation algorithm that 
produces estimates of attitude, angular rate, position, and velocity for a low earth orbit (LEO) 
spacecraft. The system relies on GPS phase, range, and range rate data as well as magnetometer 
data. The algorithm used is an extended Kalman filter (EKF), blended with a ‘pseudo-linear’ 
Kalman filter algorithm, that was developed to provide LEO attitude, orbit, and rate estimates 
using magnetometer and sun sensor data [Deut 98]. For many LEO spacecraft the sun data is 
available during only a portion of the orbit. However, GPS data is available continuously 
throughout the orbit. GPS can produce accurate orbit estimates and combining GPS and 
magnetometer data improves the attitude and rate estimates [Oshm 99]. The magnetometer based 
EKF can converge from large initial errors in position, velocity, and attitude [Deut 97], 
Combining the magnetometer and GPS data into a single EKF will provide a more robust and 
accurate system. 

The EKF is based on the existing EKF of [Deut 98]. The GPS measurement models for phase, 
range, and range rate are incorporated into the existing structure of the filter. The original EKF 
produced the orbit estimates in terms of Keplerian elements. Due to the nature of the GPS 
measurements and ease of computation, the orbit estimates are converted to the Cartesian position 
and velocity. The measurement model for the magnetometer is adjusted for this change in the 
state. 


The estimation of the angular velocity, or rate, typically involves numerical differentiation that 
introduces noise into the update of the rate estimate. An alternative approach, developed by Bar- 
Itzhack, avoids differentiation by incorporating the kinematics into the filter [Bari 00]. This 
method is included for the magnetometer only. 

2 - ALGORITHM DEVELOPMENT 

The GPS/magnetometer navigation system is based on an extended Kalman filter algorithm, 
blended with a ‘pseudo-linear’ Kalman filter [Deut 98]. The general EKF equations are: 


X(t) = f(X(t),t) + w(t) (2.1) 

Z k+1 =h k+ ,(X(t k+1 )) + n k+1 (2.2) 

The general equations for the ‘pseudo-linear’ Kalman filter are: 

X(t) = F(X(t k )) * X(t) + w(t) (2.3) 

Z|c+i = H k+] *X k+1 +n k+l (2.4) 


The output filter state is given as 


X t =[R t ,V t , C T ,q T ,o) T ,b T ] 

Where R and V are the spacecraft position and velocity vectors, respectively, C is a vector of GPS 
receiver clock errors, g is the attitude quaternion, co is the rotation rate, b is the magnetic field in 
body coordinates. 

The filter processes data from the magnetometer, a 2-baseline GPS receiver, and reaction wheels 
(used in Euler’s equations). The GPS data, as mentioned, consists of range, range-rate, and phase 
measurements. Nominally, the filter first propagates the state estimate and filter covariance and 
then performs an update with the magnetometer data followed by an update with the GPS data. 
Data from a single GPS satellite observation is included in each update. The propagation of the 
state estimate and the filter covariance will be described first. The description of the state and 
covariance update will follow, including the development of both the magnetometer and GPS 
measurement matrices. 

2.1 Filter Propagation 

The dynamics for the position and velocity, the quaternion, and the rotation rate and b are treated 
as uncoupled. The position and velocity are propagated numerically with a simple two-body 
model, using equation (2.1). The clock errors are propagated according to a two state random- 
process model [Brow 97], The quaternion is propagated numerically using the kinematic equation 
[Wert 84]. Finally, the rotation rate is propagated using Euler’s equation and b is propagated with 
the equation expressing the time rate of change of a vector with respect to a rotating coordinate 
system. The two equations are augmented as [Bari 00] 
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where: B = [b m x] = cross product matrix containing the measured magnetic field vector, bm 


I = inertia matrix 

h = angular momentum of reaction wheels 
T = external torques on the spacecraft 

U = Db r 

D = direction cosine matrix which transforms from inertial to body coordinates 

br = reference magnetic field vector in inertial coordinates, computed from the estimated 

position 

Equation (2.5) is treated as a ‘pseudo-linear’ equation in which the dynamics matrix contains the 
current best estimate of the rate, to . 


The filter covariance matrix is propagated according to 

P k+1 (-) = *„ (X k (+))P k (+)<*>/ (X k (+)) + Q k (2.6) 

■ • ■ F*At 

<f> = the state transition matrix = e 
F= df/dX 

Q = matrix of standard deviations of the white noise processes driving each of states. 

In general, the linearized error dynamics are given as 

x(t) = F*x(t) (2-7) 


Where the dynamics matrix, F, is computed for the position, velocity, and quaternion. The error 
state is represented by x for the position, velocity and quaternion. The rate and b state equation is 
linear and therefore the dynamics matrix is given in (2.5) (and is used to propagate the error in the 
actual states © and b). The matrix F is partitioned as 
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( 2 . 8 ) 


where: F orbit =5f orbit f dX 

Uit T = [Y, -p*R/|R| 3 ] 

Jl O' 

F c'ock -|o Q 

Fattitude [©x] 

F rate,b — See (2.5) 

The filter uses the ‘multiplicative’ approach in the attitude estimation [Deut 97], therefore the error 
state for the attitude is a vector of attitude errors. 


2.2 Filter Update 

The EKF state and covariance update equations are given as 


X k+1 (+) = X k+1 (-) + K k+1 [y t+i -h k+I (X k+I (-)] (2.9) 

P k+ i(+) = [I “ K k+1 H k+1 (X k+l (-)]P k+1 (-)[I -K k+1 H k+I (X k+1 (-)] T + K k+) RK k+) (2.10) 

The Kalman gain is computed according to 

K k+1 =P k+l (-)*H k+1 T (X k (-))[H k+1 (X fc (-))P k+1 (-)H k+l T (X k (-)) + R w ]- 1 (2.11) 


where H = 3h / 3X 

The measurement matrix, H, is partitioned as H = [H or bit, 1 (maudc. H ra(e ,b]. The individual 
measurement matrices are defined according to the two measurements used, namely the 
magnetometer measurement and the GPS phase, range, and range-rate measurements. First, the 
development of the magnetometer measurement matrix will be presented, followed by the 
derivation of the measurement matrix for GPS. 

Magnetometer Update 

The measurement matrix for the orbit terms, H orb i t , was developed in [Shor 95] in terms of 
spherical coordinates. The additional relationship to express the partial derivative in terms of the 
position and velocity is straightforward. Obviously, the portion of H or bit related to the clock errors 
is zero for the magnetometer. The measurement matrix for the attitude states is given in [Deut 98] 
as 


Hattitude [-b m x] 


( 2 . 12 ) 


The innovations for both the orbit states and the attitude states is given as [Deut 98] 


= b m -Db r 


(2.13) 


where D = (I - [ax])D =the direction cosine matrix computed with the filter estimated quaternion. 
D = true direction cosine matrix which transforms from inertial to body coordinates 


Finally, the measurement matrix, 1 I rat e,b> for the rate and b is defined by the following linear 
relationship [Bari 00] 



The innovations for the rate and b states is then 


— k+l — m ^1 ra te,z — rate,b,k+l 
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(2.14) 


(2.15) 


where X rate b l+1 (-) = the estimated rate and b components of the state vector, and the other 
parameters are as defined above. 


GPS Update 


The measurement matrix for the position, velocity, and clock states (combined) is a standard 
model based on range and range-rate data (actual GPS measurements would be pseudo-ranges and 
delta-ranges, but for the purposes of testing, range and range-rate measurements were generated). 
The measurement matrix is obtained by linearizing the range and range-rate equations that are 
functions of the spacecraft position and velocity, respectively. The resulting measurement 
matrices are [Brow 97] 
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(2.16) 


Where c x , c y , c z are the appropriate direction cosines between the spacecraft and the GPS satellite 
in view. The innovations are computed by differencing the measured range and range-rate with 
the computed values (computed range and range-rate are based on the filter estimated position and 
velocity), given as 
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(2.17) 


where p, p c = measured and computed range, respectively 

p,p c = measured and computed range-rate, respectively 

The measurement matrix for the attitude states is found by linearizing the relationship between the 
phase measurement and the line of sight vector to the GPS satellite in view of the antenna. 
Assuming a 2-baseline GPS receiver, the phase measurement for a given baseline can be written as 

y ji = £*jDs ri + Uj (2.18) 


i = GPS satellite designation 
j = antenna designation (j=l, 2) 
y = the phase measurement 

aj = unit vector representing the baseline direction for antenna j 
Sir = unit vector to GPS satellite i, in inertial coordinates 
Vj = white sequence measurement noise 

The computed phase measurement is given as: 


(2.19) 


The innovations, Z k +i, is computed by subtracting (2.19) from (2.18). The equation can then be 
recast to produce the necessary linear relationship between the Z k+ i and the attitude error vector, 
defining the attitude measurement matrix for the phase measurements 

z ji =y* ~ 9 » =-»![§. i.mX]a + Oi =H an a + o i (2.20) 


where s aim 


y.i 

y 2i 

V 1 - -y'i 


= measured unit vector to GPS satellite i in body coordinates 


a= attitude error state 


Finally, the measurement matrix for the rate states is developed by first differentiating the cosine 
of the angle between the baseline and the line of sight to a given GPS satellite gives (dropping the 
satellite and antenna designations) 

-sin(9)0 = a T Ds r + a T Ds r +a T Ds r (2.21) 

where the quantities are as defined above. Substituting in the kinematic equation for the direction 
cosine matrix [Wert 84], and noting that a= 0 gives 

-sin(0)0 = -a T [cox]Ds r +a T Ds r (2.22) 

Defining the ‘measurement’ to be 

y^ = a T Ds + sin(0)0 (2.23) 

and the measurement matrix, H ra t e ,b as 

H rate>b = |f)s] T [a x] 0 lx3 ] (2.24) 

results in the following linear update equation 

y ffl =H ratt , b X ra te.b + fi @.25) 

The innovations for the rate estimation are then given as 

Z = y ra -H rate b X rate b (2.26) 

where X rateb = the estimated rate and z state elements. The elements in y m , given in (2.23), are 
computed numerically. The derivative of the line of sight, s, can be computed from the given 

GPS position and velocity and the estimated position and velocity. The angle, 0, is computed 
from the phase measurement and is numerically differentiated. Alternatively, the ‘estimation 
approach’ of [Bari 00] could be utilized for the GPS rate estimation. 

3 - SIMULATION AND TEST CONDITIONS 

The test scenario consisted of a simulated orbit for a low earth orbit spacecraft and simulated orbits 
for the GPS constellation. The spacecraft rotated at a 1 revolution per orbit rate about the y body 
axis in an earth pointing attitude, at an altitude of 700 km and inclination of 98 degrees. The two 
GPS antennas are aligned along the spacecraft x and y body axes, on the anti-nadir side of the 
spacecraft. The magnetometer is aligned along the spacecraft body axes. The magnetometer and 
reaction wheel data contained an added white noise, the GPS measurements had no corruption. 
All the data was output at a rate of 1 Hz. The errors between the true state and the a priori filter 
state are given in Table I, except for the state, b. The state, b, was initialized with the first 
magnetometer measurement. The true state of the filter is output at 1 Hz also. 
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Table 1: A Priori Errors in the Filter State 


STATE 

Position 

Velocity 

Quaternion 

Rate 


ERROR 

100 km/axis 

0.5 km/sec/axis 

1 2 deg rotation error 

Orbit rate in the pitch axis only 


4 - RESULTS 

The initial test of the GP S/magnetometer navigation system included 30 minutes of simulated data. 
Figures 1 through 4 show the results of this first test. All the state parameters plotted appear to be 
converging after approximately 10 minutes. All of the states experience large deviations during 
the first 10 minutes. Figure 1 shows the position errors about each axis. By 30 minutes the errors 
are approximately 1 km in each axis. Figure 2 shows the velocity errors about each axis. Figure 3 
shows the angular attitude error about each of the body axes. By 30 minutes the errors are reduced 
to less than 1 degree in each axis. Finally, Figure 4 shows the rate errors about each axis. The rate 
errors in the x and y body axes look noisier than those in the z body axis. A plot of the additional 
state that was added according to the ‘estimation approach’ of [Bari 00] is not included. It follows 
the magnetic field vector in body coordinates, as expected. 
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Fig. 2. Velocity Errors 






u 

■o 



Fig. 3. Attitude Errors 



Fig. 4. Rate Errors 

In this preliminary test case, the filter included data from the magnetometer and only one GPS 
satellite per update cycle. Additional GPS satellites in view could easily be processed per update 
cycle. Advantages of processing only one per cycle are the reduction in processing time per 
update and in the reduction of data that must be carried between two update cycles (for use in the 
derivative computations). 

5 - CONCLUSIONS 

A combined GPS/magnetometer navigation system, based on an Extended Kalman filter algorithm, 
blended with a ‘pseudo-linear’ Kalman filter algorithm is presented. The system simultaneously 
estimates the position, velocity, attitude, and rate for a low earth orbit spacecraft using a 
magnetometer and a 2-baseline GPS antenna/receiver and reaction wheels. An initial test 
conducted with a 30 minute span of simulated magnetometer and GPS data shows that the filter 
can estimate all the state elements. Additional tests, both with simulated and eventually real 
spacecraft data, and tuning will determine the accuracy potential of this system. It is anticipated 
that this system can provide robust, accurate, and autonomous navigation estimates by relying on 
both the magnetometer and GPS. 
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